import datetime
from pathlib import Path

import matplotlib.pyplot as plt
import numpy as np
import sympy as sp
from docxtpl import DocxTemplate, InlineImage
from scipy.optimize import fsolve

# tau_relay_stabilization_1.0_formating test
print("DZ condition_1")

# численные данные
ksg = 2
ke = 11.3
Te = 0.23
a = 0.027

# начальные условия
y0 = 0.2
y10 = 0  # тут y01 = y'0

# 1. получение фазового портрета
# Интегральное уравнение фазовой траектории:
dl = -1
C1 = -5.8
x02 = sp.symbols("x02")


def phase_trajectory(x02):
    """Функция для построения фазового портрета

    ПАРАМЕТРЫ
    ---------
    dl = -1
    C1 = -5.8
    x02 = sp.symbols("x02") уравнение фазовой траектории (Интегральное)
    численные данные
    ksg = 2
    ke = 11.3
    Te = 0.23
    a = 0.027

    РЕЗУЛЬТАТЫ
    ----------
    x1 = np.arange(-10, 10, 0.01)
    x2 = [phase_trajectory(x02) for x02 in x1]
    """
    return -Te * x02 - Te * ke * dl * np.log(abs(x02 - ke * dl)) + C1


x1 = np.arange(-10, 10, 0.01)
x2 = [phase_trajectory(x02) for x02 in x1]

# вывод фазовой траектории
fig1 = plt.figure()
plt.plot(x2, x1)
plt.grid(True)
plt.show()
fig1.savefig("dz/tau_relay_stabilization/Images/image1.png")


# Линии переключения на фазовой плоскости проходят через точки:
x11 = a / ksg
print(f"Линии переключения на фазовой плоскости проходят через точки:{x11, -x11}")

# 2. Метод точечных преобразований
# уравнение, корень которого равен максимальному значению скорости:
v = sp.symbols("v")


def max_velocity(v):
    """Функция нахождения максимальной скорости

    ПАРАМЕТРЫ
    ---------
    v = sp.symbols("v")  корень данного уравнения равен максимальному значению скорости
    численные данные
    ksg = 2
    ke = 11.3
    Te = 0.23
    a = 0.027

    РЕЗУЛЬТАТЫ
    ----------
    xx = np.arange(0, 4, 0.01)
    mv = [max_velocity(v) for v in xx]
    """
    return 2 * Te * v + 2 * a / ksg + Te * ke * np.log(abs((ke - v) / (ke + v)))


xx = np.arange(0, 4, 0.01)
mv = [max_velocity(v) for v in xx]

# вывод графика
fig2 = plt.figure()
plt.plot(xx, mv)
plt.grid(True)
plt.show()
fig2.savefig("dz/tau_relay_stabilization/Images/image2.png")

# Максимальная скорость:
vmax = fsolve(max_velocity, 1.0)  # в рыбе это x2m
print(f"Максимальная скорость:{vmax}")


# Определение константы интегрирования:
def cc(x2):
    """Функция определения константы интегрирования

    ПАРАМЕТРЫ
    ---------
    ksg = 2
    ke = 11.3
    Te = 0.23
    a = 0.027

    РЕЗУЛЬТАТЫ
    ----------
    c1 = cc(vmax)
    print(f"константа интегрирования:{c1}")
    """
    return x11 + Te * x2 + Te * ke * dl * np.log(abs(x2 - ke * dl))


c1 = cc(vmax)
print(f"константа интегрирования:{c1}")


# учет константы интегрирования в уравнении фазового портрета
def phase_trajectory_f(x02):
    """Функция, которая учитывает константы интегрирования в уравнении фазового портрета

    ПАРАМЕТРЫ
    ---------
    ksg = 2
    ke = 11.3
    Te = 0.23
    a = 0.027
    """
    return -Te * x02 - Te * ke * dl * np.log(abs(x02 - ke * dl)) + c1


# траектория второго витка фазового портрета_2
def n_phase_trajectory_f(x02):
    """Функция, которая выводит траекторию второго витка фазового портрета_2

    ПАРАМЕТРЫ
    ---------
    ksg = 2
    ke = 11.3
    Te = 0.23
    a = 0.027

    РЕЗУЛЬТАТЫ
    ----------
    x1 = np.arange(-vmax, vmax, 0.01)
    x2 = [phase_trajectory_f(x02) for x02 in x1]
    nx1 = -x1
    nx2 = [n_phase_trajectory_f(x02) for x02 in x1]
    """
    return -(-Te * x02 - Te * ke * dl * np.log(abs(x02 - ke * dl)) + c1)


x1 = np.arange(-vmax, vmax, 0.01)
x2 = [phase_trajectory_f(x02) for x02 in x1]
nx1 = -x1
nx2 = [n_phase_trajectory_f(x02) for x02 in x1]

# вывод фазового портрета
fig3 = plt.figure()
plt.plot(x2, x1)
plt.plot(nx2, nx1)
plt.grid(True)
plt.show()
fig3.savefig("dz/tau_relay_stabilization/Images/image3.png")

# амплитуда автоколебаний на графике:
x1m = phase_trajectory_f(0)
print(f"амплитуда автоколебаний:{x1m}")

# определение периода и частоты автоколебаний:
T = 2 * Te * np.log(abs((ke + vmax) / (ke - vmax)))
omega1 = 2 * np.pi / T
print(f"Период:{T}")
print(f"Частота:{omega1}")

# 3. гармоническая линеаризация
# выражения для определения амплитуды и частоты автоколебаний:
AA = ((4 * ksg * ke * Te * a**2 / np.pi) ** (1 / 3)) * (1 / ksg)
omega2 = (4 * ke * a) / (np.pi * ksg * AA**2)
print(f"Амплитуда по солодовникову:{AA}")
print(f"Частота по солодовникову:{omega2}")


# 4. Графоаналитический метод
# Вычисляем вещественную (q) и мнимую (q1) части линеаризованного коэффициента передачи:


def q(A):
    """Функция, вычисляющая вещественную (q) часть линеаризованного коэффициента передачи"""  # noqa: E501
    return 4 * (1 - a**2 / (A**2)) ** (1 / 2) / (np.pi * A)


def q1(A):
    """Функция, вычисляющая мнимую (q1) часть линеаризованного коэффициента передачи"""
    return -4 * a / (np.pi * A**2)


# Уравнение гармонического баланса


def H(A):
    """Уравнение гармонического баланса

    ПАРАМЕТРЫ
    ---------
    вещественная (q)
    мнимая (q1) части линеаризованного коэффициента передачи
    """
    return (q(A) ** 2 + q1(A) ** 2) ** (1 / 2)


# ЛАфЧХ для линейного звена:


def Lwn(w):
    """ЛАФЧХ

    ПАРАМЕТРЫ
    ---------
    ksg = 2
    ke = 11.3
    Te = 0.23
    a = 0.027
    """
    return 20 * np.log(ke * ksg) - 20 * np.log(w) - 20 * np.log(Te * w)


def Fiwn(w):
    """ЛАФЧХ
    ПАРАМЕТРЫ
    ---------
    ksg = 2
    ke = 11.3
    Te = 0.23
    a = 0.027
    """
    return (-180 / np.pi) * (np.arctan(Te * w) + np.arctan(w))


# ЛАфЧХ для нелинейного звена только от амплитуды:


def LnA(A):
    """ЛАФЧХ для нелинейного звена только от амплитуды"""
    return -20 * np.log(1 / H(A))


def FiA(A):
    """Уравнение гармонического баланса

    ПАРАМЕТРЫ
    ---------
    вещественная (q)
    мнимая (q1) части линеаризованного коэффициента передачи
    """
    return -180 / np.pi + abs(np.arctan((q1(A) / q(A))) * (180 / np.pi))


# Выберем значения амплитуд для характеристики

Amplitudes_0 = [0.03, 0.05, 0.08, 0.13]


# Создаем массив значений для w

wn_values = np.linspace(-5, 100, 200)  # от -5 до 100


# Построим графики
fig4 = plt.figure(figsize=(12, 12))
# Увеличьте высоту фигуры для двух графиков

plt.subplot(2, 1, 1)  # Две строки, один столбец, первый график
plt.plot(wn_values, Lwn(wn_values))

plt.xlabel("w")
plt.ylabel("L(w)")

n_col = 0

for a in Amplitudes_0:
    color_1 = ["r", "y", "g", "b"]
    plt.axhline(y=LnA(a), color=color_1[n_col], linestyle="--")
    n_col += 1

fig4.savefig("dz/tau_relay_stabilization/Images/image4.png")

plt.show()
# Отображение графика

fig5 = plt.figure(figsize=(12, 12))

plt.subplot(2, 1, 2)  # Две строки, один столбец, второй график
plt.plot(wn_values, Fiwn(wn_values))

plt.xlabel("w")
plt.ylabel("fi(w)")

n_col = 0
for b in Amplitudes_0:
    color_1 = ["r", "y", "g", "b"]
    plt.axhline(y=FiA(b * ksg), color=color_1[n_col], linestyle="--")
    n_col += 1


plt.subplots_adjust(hspace=2)  # Установите желаемое расстояние между графиками
plt.tight_layout()
plt.show()

fig5.savefig("dz/tau_relay_stabilization/Images/image5.png")


# Создание папок, если они не существуют
images_folder = "dz/tau_relay_stabilization/Images"
Path(images_folder).mkdir(parents=True, exist_ok=True)

output_folder = "dz/tau_relay_stabilization/Otchet_docx"
Path(output_folder).mkdir(parents=True, exist_ok=True)

# Создание папки для docx файла
Path(output_folder).mkdir(parents=True, exist_ok=True)


# Далее будет оформление docx файлов.


doc = DocxTemplate("dz/tau_relay_stabilization/Otchet_docx/my_word_template.docx")
context = {
    "title": "Исследование релейной стабилизации угла крена",
    "year": datetime.datetime.now().strftime("%Y"),
    "image1": InlineImage(doc, "dz/tau_relay_stabilization/Images/image1.png"),
    "image2": InlineImage(doc, "dz/tau_relay_stabilization/Images/image2.png"),
    "image3": InlineImage(doc, "dz/tau_relay_stabilization/Images/image3.png"),
    "image4": InlineImage(doc, "dz/tau_relay_stabilization/Images/image4.png"),
    "image5": InlineImage(doc, "dz/tau_relay_stabilization/Images/image5.png"),
    "text": "Москва",
}

doc.render(context)
doc.save("dz/tau_relay_stabilization/Otchet_docx/Отчет_питон.docx")
